Advancing thermal performance through vortex generators morphing

The design of rigid vortex generators (RVG) influences the thermal performance of various technologies. We employed Discrete Adjoint-Based Optimization to show the optimal development of vortex generators. Under turbulent flow conditions, different bi-objective functions on the RVG design were examined. Specifically, we aimed at an optimal RVG shape that minimizes the pressure drop and maximizes the local heat transfer in a rectangular channel. We show that an optimal design of an RVG can be obtained using computational fluid dynamics in conjunction with the Pareto Front at a computational cost of the order ~\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$O(10^{-1})$$\end{document}O(10-1). We obtained three essential vortex generator shapes based on the RVG morphing technique. Compared to the baseline geometry of a delta winglet pair DWP, the first morphed design reduced the pressure drop by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$39\%$$\end{document}39%, however, at the expense of a \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$21\%$$\end{document}21% reduction in the Nusselt number. The second vortex generator design enhanced the heat transfer by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$18\%$$\end{document}18%, however, at the cost of a significant increase in pressure drop of about \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$40\%$$\end{document}40%. The final morphed design achieved the highest thermal performance factor of 1.28, representing a heat transfer enhancement of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$6\%$$\end{document}6% with a moderate increase in pressure drop of about \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$13\%$$\end{document}13% compared to DWP vortex generators. Furthermore, we investigated the effect of introducing different size holes on the mass reduction of vortex generators and their thermal performances. The mass of vortex generators can be reduced by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$9\%$$\end{document}9% and with an increase of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$7\%$$\end{document}7% in thermal performance factor concerning the DWP baseline. The findings of this study will lead to highly efficient lightweight heat exchangers.

Several researchers in the past performed computational modelling of vortical structures. For example, Danaila and Hélie 11 investigated the post-formation evolution of a laminar vortex ring. Sun et al. 12 studied the wake structure of a micro-ramp vortex generator in hypersonic flows. Tian et al. 13 investigated different types of longitudinal RVGs in a flat plate channel flow to increase the heat transfer through the Nusselt number and reduce the pumping power through the friction factor. Moreover, researchers conducted numerical investigations of RVG designs like annular winglet; wave-element RVG 14,15 ; curved winglet RVG with or without local perforations 16,17 ; twisted and helical RVGs 18,19 ; wavy-fins and oval-tube-bank with mounted RVGs 20,21 and longitudinal vorticity RVGs 22 . Zaman et al. 23 tried experimentally controlling an axisymmetric jet by using vortex generators. Behfard and Sohankar 24 characterized the local flow structures for a fluid flow in the presence of delta winglet pair (DWP) in a finned circular heat exchanger under a turbulent regime with a steady SST k − ε turbulence model. Oneissi et al. 25,26 investigated the influence of a new RVG design as an inclined-projected-winglet-pair (IPWP) on the heat transfer enhancement under a turbulent flow regime. Finally, Hamidouche et al. 27 conducted flow characteristics downstream delta-winglet vortex generator experiments using the stereoscopic particle image velocimetry (PIV) technique in a rectangular channel but without heat transfer. Experimental measurements of RVGs with coupled fluid flow and heat transfer are scarce in the literature, and the last works go back to the 1990s by Fiebig and Tiggelbeck 28,29 . The difficulty is to cover a wide range of operating conditions, such as different Reynolds and Prandtl numbers, RVG designs and angles of attack.
The optimization of thermofluid components by employing advanced optimization techniques in CFD, e.g., heat transfer enhancement, pressure drop, and overall weight/mass reduction, has emerged as a topic of research and development. For example, Kim and Choi 30 realized shape optimization of dimples RVGs in a rectangular channel by employing a Surface Response Optimization (SRO) method. Song et al. 31 developed an advanced optimization technique in CFD using a Surrogate-Based Algorithm (SBA) for designing a heat sink. Recently, Karkaba et al. 3 presented multi-objective function optimization through a large space exploration design developed to find optimal RVG design. They defined seven different design variables to draw/find the shape of the optimal RVG. A new enhanced RVG design was found recently by 3 through parametric shape optimization. The optimal design 3 led to a substantial increase in the thermal enhancement factor, though at the expense of computational costs. Moreover, Oh and Kim 32 studied the thermal and flow characteristics of the rectangularwinglet (RW); delta-winglet-upstream (DWU); and delta winglet-downstream (DWD) curved vortex generators (CVGs). These designs varied regarding the position angle and the radial distance in the fin-tube heat exchanger configuration. Most of their results suggested that conventional VGs are superior in performance compared to curved vortex generator CVGs. Carpio and Valencia 33 investigated the effect of introducing longitudinal vortex generators (LVGs) in arrays consisting of 9 alternating, 18 alternating, 18 aligned, 27 alternating and 39 alternating LVGs. Results show that the design with the best thermal performance corresponds to geometry with 39 LVG, with a 52% increase compared to the flat fin geometry. The case with an enhanced performance corresponds to the geometry with 18 vortex generators, which presents a performance increase of 38% to the baseline design. Therefore, it performs similarly to the 39 LVG case, despite having only half the vortex generators. Khan et al. 34 adopted artificial neural networks (ANNs) for predicting the efficiency of a double-pipe heat exchanger which employs T-W tape inserts with different wing-width ratios. The ANNs predicted with high accuracy thermal parameters such as the friction factor, the Nusselt number and the thermal performance factor (TPF) of a double pipe heat exchanger. More recently, Khan et al. 35 used ANN to optimise the performance of a double pipe heat exchanger, i.e., to obtain better heat transfer enhancement with a low-pressure drop penalty. They found a maximum value of TPF equal to 1.44 by varying five input parameters related to wing-width ratio, pitch ratio, attack angle, Reynolds number, and tube length (L).
As seen in the research mentioned above, the design method for vortex generators typically revolves around a set of geometrical parameters such as the angle of attack, inclination angle, span, height, and different shapes such as delta wing and rectangular wing. Since the parametric space is ample, it becomes a tedious task to define a concrete parametric optimization problem without considering that designing vortex generators within these parameters is somewhat restrictive. Optimization processes typically require extensive CFD calculations to explore the variations of just a few parameters. Therefore, effective optimization techniques are needed to rapidly advance the development of vortex generators' performance. To this end, gradient-based optimization methods aim to find 'sensitivity gradients' from the solution of the flow field. These will indicate how the performance will change by modifying the vortex generator's geometry. Such advanced techniques, therefore, minimize the computational cost of the optimization process. Gradient-based methods have been developed mainly for aerospace applications to enhance lift and reduce drag for airfoils geometries 36,37 but also have great potential for heat exchangers technology. The present work chooses the Adjoint method for several reasons. To the author's knowledge, Adjoints have not been used in conjunction with vortex generators in current literature. Still, they can be powerful and efficient tools in developing vortex generator-based heat exchangers. Furthermore, popular CFD codes such as ANSYS Fluent have an Adjoint solver already implemented so that the methods developed here are widely accessible to the general CFD/Thermal engineering community working on vortex generators. This will increase the impact of this research and shed more light on the Adjoint method's applicability to heat enhancement techniques based on the vortex generator technology.
In this study, we employ an adjoint-based shape optimization method in conjunction with CFD for the optimal and rapid design of RVGs in a turbulent flow regime. We define a bi-objective Adjoint-Based Optimization problem to compute the sensitivity maps by: (1) locally morphing the mesh and deforming the RVG surface topology in a controlled volume sub-domain; (2) creating the overall optimization algorithm to act in an enhanced and more rapid optimal line-search direction using the steepest descent algorithm, thus reducing the overall computational time of the Pareto Front by an order ~ O(10 −1 ) . For the first time, we show and quantify www.nature.com/scientificreports/ the direct influence of the mathematical formulation of the BOF on the Pareto Front, the overall computational time and the RVGs design optimality.

Results
The Pareto Front for the RVG designs fits well with a Pareto function (Fig. 1). The solid squares in Fig. 1 represent the adjoint-based designs obtained with a weighted sum of function (WSM) f 2 (Eq. 5) while the solid circles represent the adjoint-based RVG designs obtained with thermal performance factor (TPF) morph function f 1 (Eqs. 4,9). The design space using f 2 is much larger than the design space using f 1 represented by the coloured elliptical schematic zones in Fig. 1. Therefore, using the ABO algorithm with BOF in the form of f 1 (TPF morph) can find optimal designs more efficiently and rapidly than BOF in the form of f 2 (WSM Morph). The overall computational time to obtain the TPF morph-based optimum ( TPF = 1.28 ) design is reduced tenfold compared to the WSM morph-based optimum design ( TPF = 1.26 ). Increasing the TPF above one, even by a small percentage of 1% , gives a major improvement when optimizing a single RVG in a channel flow. A minimal increase in the TPF value, e.g., 1% , can have a much larger enhancement in a compact heat exchanger (HEX) overall performance due to thousands of locally distributed RVGs on the HEX channel surfaces. The RVG's shape evolution during the ABO process (algorithm of Fig. 10) is shown in Fig. 2. We also show the frontal area evolution in mm 2 and the different RVG designs initially at a it = 0 , after the first adjoint iteration a it = 1 and at the end of the overall iterations a it = n . From Fig. 2, we observe that the ABO algorithm, using TPF morph BOF f 1 , succeeds in deforming the RVG surface in a way that promotes the local creation or generation of additional vortices. This can be seen from the sharp edge at the upper lateral corners of the optimal RVG. Moreover, the frontal area is reduced from 220 to 197 mm 2 in optimal designs obtained by f 1 and f 2 compared to the initial DWP design.
Here we give specific attention to optimizing the mass of the optimal RVG design. As shown in Fig. 3, we investigate the mass reduction of the new RVG design by introducing different holes at a time, each of diameter d hole . Compared to DWP RVG, it is found that we can reduce the mass as in the new perforated RVG by 9.5% with a TPF increase of about 7% , only at the scale of a single RVG (see PWP-ALI-3 in Fig. 3). Nevertheless, more significant mass reduction is associated with a lower increase in the TPF, which is seen in Fig. 3, where the mass can be reduced by up to 17.5% , however, with a smaller increase of 3.5% in TPF when compared to the DWP RVG. So, more HEX designs can be produced at a larger scale using the new PWP-ALI RVG design within a turbulent flow in rectangular channels. As a result, the future HEX can have highly reduced overall mass and a good increase in its overall TPF. This is a significant finding in research and development of turbulent fluid flow in rectangular channels and its effect on designing optimal compact heat exchangers/reactors of reduced  4,9); f 2 (WSM Morph based, Eq. 5). The overall computational time to obtain the TPF morph-based optimum design is reduced tenfold compared to the WSM. This reduction in computational time can be explained by the reduction of the design space of the exploration represented schematically by the two coloured zones in red (for f 1 ) and blue (for f 2 ). The optimal RVG design obtained by f 1 has TPF = 1.28 and the optimal RVG design obtained by f 2 has TPF = 1.26. www.nature.com/scientificreports/  www.nature.com/scientificreports/ mass versus increased overall TPF. By taking a plate-fin heat exchanger with inserted vortex generators made of Aluminum material that has a density of 2700 kg/m 3 , the mass of one pair of DWP RVGs will be equivalent to approximately 2.16 g. Based on the results of Tiggelbeck et al. 38 where two rows of DWP RVGs can be inserted along a channel length of 15H for effective heat transfer enhancement. As a result, this will correspond to a mass per volume ratio of 7200 g/m 3 , i.e. 7200 g of DWP RVGs mass per one cubic meter of heat exchanger volume. The results for the amount of mass reduced per unit cubic meter of heat exchanger volume for the three perforated designs PWP-ALI-1,2 and 3 are summarised in Table 1. Figure 3b-d illustrate the local surface deformation induced by the adjoint-based optimization algorithm to maximize the TPF and minimize the overall mass. One can observe that the initial geometry (DWP) is locally deformed at different height positions from y = 1 mm to y = 19 mm concerning the bottom wall at y = 0 where the RVG is fixed. Between y = 1 mm and y = 10 mm , which corresponds to half the channel's height, the ABO algorithm deforms the surface of the DW RVG towards the lateral side away from the centre-line of the main flow. Above y = 12 mm and close to the upper wall, the ABO algorithm deforms the initial RVG design towards the centre-line of the main flow. This reversed curvature of the RVG surface has an essential feature in fluid flow physics. It permits the production of additional local vortices, thus enhancing heat recovery from the lateral walls at a reduced pressure drop in the channel. Flow structure and temperature distribution are also analysed to understand better the heat transfer enhancement mechanism by local RVGs. Figure 4 shows the velocity streamlines along 15 positions from P1 to P15. The streamlines are computed in 15 different xy-planes distributed with equal spacing between the inlet ( z = 0 , P1) and outlet ( z = 15h , P15). Between P1 and P7 in all the RVG designs, two major vortices are created downstream of the RVG. Looking at the plane P7, the two vortices generated by the BOF f 1 are clearly wider than those produced by the DWP RVG and the two extreme conditions with f 2 ( a = 1; b = 0 ) and f 2 ( a = 0; b = 1 ). Between P9 and P15, the optimal RVG design obtained with f 1 develops six local longitudinal vortices. The latter have vertical centres closer to the channel's lateral boundaries; see the red-coloured solid circles in Fig. 5, which represent the local dynamics of the vortical centres from inlet to outlet. Compared to the DWP design, the optimal RVG obtained with f 1 produces three major vortical structures more displaced downstream in the XY lateral directions. This removes more heat from the walls by forced convection and at a reduced pressure drop.

Discussion
To better investigate the turbulence production rate in the channel due to the RVGs, we compute in Fig. 6 the 2 -criterion. At P15 with optimal f 1 , it can be observed that six major vortices have better conserved strength and more uniform spanwise distribution compared to DWP and the two extreme conditions when using the BOF f 2 with a = 1 and b = 0 and with a = 0 and b = 1.
In terms of quantitative analysis of fluid flow physics and heat transfer, we compute the TPF (see Eqs. 4,9), global Nusselt number Nu g and global friction factor f g . Figure 7 clearly shows that the TPF value is increased by about 8.6% when using a WSM morph-based BOF f 2 , while the TPF value is increased by 10.3% when employing a TPF Morph based BOF f 1 . Moreover, the TPF Morph optimal design is obtained at a computational cost reduced by a factor of 10, compared to the WSM optimal design. This illuminates the advantages of employing ABO with BOF of the form of f 1 (Eqs. 4,9) rather than of the form of f 2 (Eq. 5).
For further local analysis of the flow in the channel, the readers may refer to the supplementary information document attached to the manuscript.
In conclusion, we employed the Discrete ABO method coupled with advanced CFD to optimise a DWP RVG under a turbulent flow regime at Re ≈ 4600 . The objective is to increase the heat transfer and reduce the pressure drop in a rectangular channel. In other words, we seek to improve the TPF considerably. For that, we formulated two forms of a BOF with attention to the overall mass of the RVG: (1) a TPF Morph function; (2) a weighted summation method or function (WSM) with two weighting factors. We quantitatively showed how one could obtain an optimal design of a local RVG in a channel on a Pareto Front with a very well-reduced computational time of order ~ O(10 −1 ) . The above is achieved using a TPF morph BOF within an ABO method in CFD. We showed how the ABO algorithm, coupled with the CFD solver, conducted a local surface morphing of the initial RVG design with more local surface deformation at the corner side of the VG. This ensured that the new RVG produced more local vortices towards the lateral sides of the channel. As a result, the ABO algorithm converged toward a new design of RVG, which increased the thermal performance factor (TPF) by 10.3% . Furthermore, we investigated the mass reduction of this new RVG by adding a hole. Different hole diameters were investigated. Compared to a DWP RVG, we found that the mass of the new Perforated Winglet Pair RVG (named "PWP-ALI") is reduced by 9.5% , with TPF increased by about 7% . The above implies that, at larger scales, new HEX designs can be produced using PWP-ALI RVG so that the overall HEX can have a primarily reduced mass (presence of thousands of RVGs). An essential increase in thermal performance will accompany this compared to the DWP www.nature.com/scientificreports/  4,9), Optimum design using f 2 with a = 1 and b = 0 (Eq. 5) and Optimum design using f 2 with a = 0 and b = 1 (Eq. 5).

Figure 5.
Local dynamics of the vortical centers from inlet to outlet for the initial DWP design and the optimized design obtained using a TPF Morph objective function f 1 (see Eqs. 4,9). The colour intensity of the vortical centres decreases progressively as going in the stream directly from the inlet to the outlet.   www.nature.com/scientificreports/ sen because it corresponds to a maximum TPF value for the RVG compared to other higher or lower values of Re, as was shown in the literature by Oneissi et al. 25 . The conjugated fluid dynamics and heat transfer are solved through the following consecutive mass, momentum and energy conservation equations: where U is the velocity vector, T the temperature, ρ the fluid density, µ the fluid's dynamic viscosity, µ t the turbulent dynamic viscosity, Cp the specific heat, Pr and Pr t are the Prandtl and turbulent Prandtl numbers, respectively. The operating fluid is air. The Reynolds averaged Navier-Stokes (RANS) in conjunction with the k − ω − SST (Shear-Stress transport) model 39 are employed to solve for the turbulent flow field. The SST model includes two additional partial differential equations coupled to 1, 2 and 3 to account for the turbulent kinetic energy k and the turbulent dissipation rate ω . The SST model employs a k − ω model formulation in the inner parts of the boundary layer and the k − ε model in the freestream. The k − ω − SST model has shown accurate predictions and numerical convergence for flows around RVGs 26 .
The boundary conditions ( Fig. 9) are set as follows • Inlet: A constant temperature T in and uniform velocity at the inlet U in .
• Upper and lower channel surfaces: a constant temperature T w .
• Left and right channel surfaces: a symmetry boundary condition.
• VG surface: adiabatic wall with zero heat flux.
• Outlet: a zero-pressure Neumann boundary condition.

Adjoint-based optimization in CFD. The flow chart of the adjoint-based optimization (ABO) algorithm
is shown in Fig. 10, and the components of ABO are discussed below.
Discrete adjoint method. Adjoint optimization methods require no parameterization of the geometry, which is the case in other methods, such as in parametric optimization techniques in CFD 3 . Moreover, adjoint-based methods can produce better engineering designs at an overall reduced cost. Compared to the continuous adjoint www.nature.com/scientificreports/ method, the discrete approach computes the sensitivity field more accurately, especially in turbulent flow conditions.
Optimization problem. Two mathematical problems of adjoint-based optimization are formulated and solved by minimizing two different global bi-objective functions (BOF) f 1 and f 2 as the following: where f g is the global friction factor, Nu g the global Nusselt number, and the subscript "0" denotes an empty channel. f 1 is a linearly weighted sum morph (WSM) objective function with two real parameters a and b such that a, b ∈ [0 − 1] . f 2 is a ratio form single objective defined as a TPF Morph based function. The global Nusselt number is calculated using the logarithmic mean temperature difference method as follows: where q ′′ w is the average wall heat flux at the channel top and bottom walls, D h is the hydraulic diameter of the channel and is equal to 2h and k is the thermal conductivity of the flowing fluid.
The log mean temperature difference LMTD is calculated using the following equation: www.nature.com/scientificreports/ where T i is the temperature difference between the wall and the average temperature at the channel's inlet and T o is the temperature difference between the wall and the average temperature at the channel's outlet, respectively. Moreover, the global friction factor f g is related to the pressure drop p between the inlet and the outlet of the channel by: where U in is the air velocity at the inlet of the channel. The thermal peformance factor (TPF) of Eq. (4) is given by: To analyze the optimized RVGs performance, the local Nusselt number inside the channel is computed as a function of the dimensionless local position z/h in the flow direction: where and where T z is the local temperature, T w is the wall temperature, D h the hydraulic diameter, q ′′ z is the local wall heat flux, Nu z the local Nusselt number and h z the local heat transfer coefficient.
The local friction factor (or fanning friction factor) in the channel f z is computed as the following: where τ w is the shear stress at the wall. Finally, the Reynolds number is calculated based on the hydraulic diameter of the channel as follows: Sensitivity maps and mesh morphing control. For computing the shape sensitivity of the RVG concerning the objective function F such that F ∈ [f 1 , f 2 ] , we denote δ as the Cartesian coordinates (X, Y, Z) locations for every cell or node in our computational domain in the control volume v = I A × w A × h (see Fig. 8) enveloping the RVG surface. Thus, δ is an input vector to the overall adjoint optimization problem such that: The first term dF dδ represents the total sensitivity of F with respect to δ at a given cell center or node. This is required to perform the mesh morphing and apply an optimal iterative displacement to the RVG surface (in outward or inward direction, depending on the sign) to satisfy the minimization of F. A mesh morphing algorithm guarantees the quality of the morphed mesh through local control of the positive volume of each finite volume cell. The steepest descent algorithm guarantees the convergence of the optimization problem toward a local minimum of F in the design space. The second term ∂F ∂δ corresponds to the change in F due to a change in the position at a given cell center. The final term ∂R ∂δ is the change in F due to the sensitivity of the flow solution to changes in the www.nature.com/scientificreports/ local position of the RVG. R is also known as the Residuals vector/matrix of the flow solution, and T m is a vector of Lagrange multipliers, representing the adjoint solution variables, that is found by satisfying: Equation 16 is also known as the adjoint problem that is solved using the same fluid flow solvers, and β is the vector of fluid flow solution variables.
Polyhedral cells are generated in the computational domain employing a local mesh refinement method and a volume control imposed in the VG design space zone. The automatic local refinement algorithm ensures that the value of the dimensionless wall distance y + in the wall-bounded channel flow always satisfies y + ≤ 1 during the iterative procedure of mesh morphing.
Several polyhedral grids have been generated to choose the appropriate mesh size in a mesh sensitivity analysis procedure. The grid-convergence index (GCI) chooses the appropriate initial mesh size. A polyhedral mesh which contains 1.3 × 10 6 cells is adopted with less than 2% GCI value 40 validated with the experimental DWP design of Tiggelbeck et al. 29 at same Reynolds number Re = 4600 . The comparison between numerical simulations and experimental results is summarized in Table 2 based on the calculation of the Nusselt number normalized by the Nusselt number corresponding to an empty channel and a normalized friction factor and the thermal performance factor (TPF). The results in the table show that the numerical simulation model agrees with the experimental results done by Tiggelbeck et al. 29 with an error that is less than 3.5%.
Another approach to verify the choice of a numerical model and mesh size is the validation with fully developed turbulent flow in an empty channel. The numerical results are compared with the correlation of Dittus-Boelter 41 for the evaluation of the global Nusselt number and the correlation of Blasius 42 for the assessment of the global friction factor. The CFD results for an empty channel compared to the empirical correlations show an error of less than 2% for both the global Nusselt number and the friction factor.
In each ABO inner loop, the convergence of the CFD solvers is verified by looking at the friction factor f and Nusselt number (Nu) stability versus sub-iterations.
A mesh strategy in each ABO outer loop is set as follows: a control volume v is defined as an internal subdomain (Fig. 8) for the local deformation (mesh morphing) of the RVG surface based on the iteratively computed sensitivity map/field dF dδ (see Fig. 10).

Data availability
Data is available upon request from the corresponding authors.